We are IntechOpen, the world's leading publisher of Open Access books Built by scientists, for scientists

Open access books available 5,300

130,000 155M

International authors and editors

Downloads

Our authors are among the

most cited scientists 154 TOP 1%

Selection of our books indexed in the Book Citation Index in Web of Science™ Core Collection (BKCI)

# Interested in publishing with us? Contact book.department@intechopen.com

Numbers displayed above are based on latest data collected. For more information visit www.intechopen.com

# **Topological Characterization and Advanced Noise-Filtering Techniques for Phase Unwrapping of Interferometric Data Stacks**

# Pasquale Imperatore and Antonio Pepe

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/61847

### **Abstract**

This chapter addresses the problem of phase unwrapping interferometric data stacks, ob‐ tained by multiple SAR acquisitions over the same area on the ground, with a twofold objective. First, a rigorous gradient-based formulation for the multichannel phase un‐ wrapping (MCh-PhU) problem is systematically established, thus capturing the intrinsic topological character of the problem. The presented mathematical formulation is consis‐ tent with the theoretical foundation of the *discrete calculus*. Then within the considered theoretical framework, we formally describe an innovative procedure for the noise filter‐ ing of time-redundant multichannel multilook interferograms. The strategy underlying the adopted multichannel noise filtering (MCh-NF) procedure arises from the key obser‐ vation that multilook interferograms are not fully time consistent due to multilook opera‐ tions independently applied on each single interferogram. Accordingly, the presented MCh-NF procedure suitably exploits the temporal mutual relationships of the interfero‐ grams. Finally, we present some experimental results on real data and show the effective‐ ness of our approach applied within the well-known small baseline subset (SBAS) processing chain, thus finally retrieving the relevant Earth's surface deformation time ser‐ ies for geospatial phenomena analysis and understanding.

**Keywords:** SAR interferometry, phase unwrapping, discrete calculus

## **1. Introduction**

Multichannel (or multitemporal) InSAR techniques address the processing of interferometric data stacks obtained by combining multiple SAR acquisitions over the same area. These approaches can be essentially categorized in two main classes: persistent scatterers (PS) and

© 2016 The Author(s). Licensee InTech. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

small baseline (SB)-based techniques. The solution of the multichannel phase unwrapping (MCh-PhU) problem is generally required in multichannel InSAR techniques, whenever multidimensional SAR data set, conveying information about complex Earth's crust events, have to be systematically investigated on suitable space-time scales for geospatial phenomena understanding [1–29]. In this chapter, we focus on two different related main issues.

Primarily, we present a rigorous gradient-based formulation of the MCh-PhU problem that is consistent with the theoretical foundation of the discrete calculus [30–34]. Emphasis is placed on the topological characterization of the underlying discrete setting provided by the *differen‐ tial operators* of the discrete calculus, which are formally associated with matrix operators and represent the discrete counterparts of the classical differential operators of the vector calculus. Accordingly, MCh-PhU problem formulation is systematically established in terms of discrete differential operators, which are defined by the topology of the intrinsically discrete spaces upon which they act, thus capturing the essential topological character of the problem within a systematic matrix formalism [35]. It is worth highlighting that our approach provides an unambiguous and theoretical-consistent formalism for the MCh-PhU problem, overcoming the conceptual inconsistencies of the existing gradient-based formulations [1, 17, 29]. Indeed, the existing approaches pose some conceptual limitations from a mathematical viewpoint since they rely on an intrinsically discrete setting based description and, at the same time, resort to the concepts of gradient and curl of the conventional vectorial calculus, which inherently imply a reference to an underlying continuum space and the notion of the infinitesimal [30]. In addition, the proposed formal framework enables meaningful analytical investigations on a mathematical consistent playground, also providing interesting implications and permitting to express previous obtained results in a more general way.

Then we present an innovative procedure to filter out the noise affecting the phase components of a redundant set of (multitemporal) multilook small-baseline interferograms. This is achieved by independently solving, for each pixel of the scene, a *nonlinear optimization* problem based on computing the wrapped phase vector that minimizes the (weighted) circular variance of the difference between the original and noise-filtered interferograms [43]. This noisefiltering procedure arises from the key observation that multilook interferograms are not fully time consistent because they are generated through multilook operations that are independ‐ ently carried out on each single interferogram. Indeed, the wrapped discrete curl of the interferometric phases defined on a graph whose nodes and edges describe SAR acquisitions (in the time/perpendicular baseline domain) and inherent interferograms, respectively, is different from zero. This modulo-2π cyclic inconsistency of multichannel interferometric phases is properly handled by the presented multichannel noise-filtering (MCh-NF) proce‐ dure. The presented technique is very easy to implement because it does not imply a prelimi‐ nary time-consuming selection of statistically homogenous pixels (SHP), as for instance required by the *SqueeSAR* technique [44], and it has no need of any *a priori* information on the statistics of complex-valued SAR images. The effectiveness of the presented noise-filtering approach as well as its impact on the quality of multichannel phase unwrapping procedures are also fully investigated.

# **2. Multichannel phase unwrapping problem**

In this Section, we review the mathematical formulation of the multichannel phase unwrap‐ ping (MCh-PhU) problem within the purview of the discrete calculus. As a matter of fact, a graph-based description naturally arises in formulating the MCh-PhU problem due to the underlying discrete irregular data structure. Indeed, as far as discrete settings (e.g., graphs) are concerned, resorting to the conventional vectorial calculus might not be adequate since it inherently implies a reference to an underlying continuum space. On the contrary, *discrete calculus* offers a rigorous methodological framework since it treats a discrete domain as entirely its own entity. In particular, discrete calculus provides proper differential operators that make it possible to purely operate onto a finite, discrete structure without referring to the continuous space and notion of the infinitesimal [30]. More specifically, the introduction of some wellknown algebraic structures [30–33] capturing the essential topological character of the underlying graphs permits to phrase pertinent *differential operators* as matrices. Therefore, one of the most important consequences of this approach is that the purely topological nature of the discrete differential operators is made more apparent and concrete. Accordingly, by systematically adopting the relevant key concepts and tools available within this theoretical framework, we here provide a description of the MCh-PhU problem on a suitable mathemat‐ ical playground. For such a purpose, we first establish the notation and terminology used throughout the subsequent Sections 2.1, and then the problem at hand is reviewed within this formalism (Sections 2.2 and 2.3). We remark that the focus here is on presenting key concepts that are useful for the following analyses; however, a comprehensive treatment of the discrete calculus and related huge fields of mathematics (e.g., algebraic topology, exterior calculus, and differential forms) is clearly beyond the scope of this work but can be found in refs. [30–33].

## **2.1. The theoretical framework of discrete calculus**

A graph G(V ;ℰ) is defined by two sets: V and ℰ. The former is the set of nodes (or vertices) of the graph, and the latter represents the corresponding set of edges. Let *Q* and *M* be the cardinality of V and ℰ, respectively. The vector space ℝ *<sup>M</sup>* is referred to as the *edge space*, and the vector space ℝ *Q* is referred to as the *vertex space*, with ℝ denoting the field of real numbers. Without loss of generality, we here assume that the graph G is *connected* (i.e., every pair of vertices in the graph is connected [33]). Moreover, an orientation establishes a default direction on an edge that is considered positive or negative, thus yielding an oriented graph. The *M* × *Q incidence matrix Π* = *Πmq* of an oriented graph G specifies its edge–node connectivity relations, whose entries are defined as follows [30–33]:

$$
\Pi\_{mq} \equiv \begin{cases} -1, & \text{if } m \text{th edge starts at } q \text{th node} \\ +1, & \text{if } m \text{th edge ends at } q \text{th node} \\ 0, & \text{otherwise} \end{cases} \tag{1}
$$

with *m* = 1, 2,..., *M* and *q* = 1, 2,..., *Q*. It is important to note that the column rank of *Π* is *Q* − 1. The incidence matrix *Π* generates an orthogonal decomposition ℝ(*Π*) ⊕ N (*Π* T )=ℝ*<sup>M</sup>* , where ℝ(*Π*) is the column space of the incidence matrix *Π*, and N (*Π* T ) denotes the kernel (or nullspace) of the matrix *Π* T . The notion of cycle space is also fundamental in graph theory. The *cycle space* of the graph G, namely, C =C(G), is the subspace of edge space ℝ *<sup>M</sup>* spanned by all the cycles in G. The dimension of C(G) is also referred to as the *cyclomatic number* of G [33]. It is also well known that, for every *connected* graph G with *Q* nodes and *M* edges, the dimension of the cycle space is given by *R* = dim(C (G)) = *M* − *Q* + 1 [30]. Each basis for C(G) (i.e., the cycle basis) is therefore uniquely specified by an *M* × *R* matrix *Ω*, called *cycle matrix*. Thus, the column vectors of *Ω* = **ω** 1 , …, **ω** *R* form a basis for an *R*-dimensional vector subspace (the cycle space of C(G) of ℝ *<sup>M</sup>* . C(G) is indeed the null-space of *Π* T so that a cycle basis provides a basis for N (*Π* T ) [30]. Accordingly, a fundamental property of a linear graph is expressed by the remarkable relations:

$$
\Pi^\top \mathfrak{Q} = \mathfrak{0} \tag{2}
$$

$$\mathbf{Q}^{\mathsf{T}} \amalg = \mathbf{0} \tag{3}$$

Indeed, several methods for defining a cycle set have been studied, and they can be used to define incidence relations between edges and cycles. Specifically, the definition of a cycle set from the edge set can be obtained algebraically and geometrically (i.e., from an embedding). Algebraic methods find a suitable *M* × *R* matrix *Ω* whose columns provide a basis for the nullspace of *Π* T , with *R* =dim(N (*Π* T ) ). Geometric methods for defining a cycle set (i.e., from an embedding) permit to identify algorithmically a cycle set (representing the faces) in this embedding and may be used to produce the edge–face incidence matrix *Ω* (as illustrated in Figure 1). In particular, it is possible to consider a normalized irreducible cycle basis forming elementary (or irreducible) cycles [30–33], i.e., cycles that contain no other cycles, so that we can associate to each elementary cycle an elementary cycle vector **ω** *r* = *ω*<sup>1</sup> , ..., *ω<sup>M</sup>* T , whose entries are defined as follows:

Accordingly, the so defined *Ω* = **ω** 1 , …, **ω** *R* provides a particularly attractive basis for N (*Π* T ), i.e., the cycle basis formed by all the elementary cycle vectors associated with the elementary cycles in G. We also note that *Ω* defines the incidence connectivity relations between edges and cycles (see Figure 1).

It is instructive to highlight that the topological operators *Π*, *Π* T , and *Ω* provide the discrete counterparts of the classical *gradient* (∇), *divergence* (∇ ⋅ ), and *curl* (∇ × ) operators of the vector calculus for continuous space, respectively. Accordingly, they can be regarded as differential operators on the discrete setting [30]. In addition, it is worth emphasizing that identities (2) and (3) mimic the properties of their classical vector calculus counterparts ∇ ⋅∇ × =0 (*div curl* = 0) and ∇ ×∇ =0 (*curl grad* = 0), respectively. It should be pointed out that *Π* yields differences along edges of nodal "potentials" represented by *Π***x**. Conversely, given an arbitrary **f**∈ℝ*<sup>M</sup>* , a solution of the equation *Π***x**=**f** (if it exists) is called the *potential* of **f**. Note also that **x** (if it exists) is not unique since the constant column <sup>=</sup> 1, ..., 1 <sup>T</sup>∈ℝ*<sup>Q</sup>* is an element of the kernel of *Π*. Of course, not every **f**∈ℝ*<sup>M</sup>* is the discrete gradient of some **x** since **f** may contain a curl component. Indeed, a prescribed **f**∈ℝ*<sup>M</sup>* can be written as a nodal difference (**f**=*Π***x**) if it is *cyclically consistent*, i.e., if it satisfies *Ω* <sup>T</sup> **f**=**0** (i.e., there is no component of the flow in the cycle space). Note also that *Ω* is the (cross) differential operator of the graph whose expression can be given in terms of a normalized cycle basis; N (*Ω* <sup>T</sup> ) denotes the subspace of ℝ *<sup>M</sup>* with zero flow circulation (curl-free) around cycles. Moreover, *Π* T **f** yields nodal accumu‐ lations from flows along edges. As a result, the differential operators, as basic tools of the discrete calculus, have been established and properly phrased on the discrete space. This mathematical abstraction meaningfully captures the topological structure of the underlying discrete setting. Note that the topological characterization of the graph is essentially embodied in the algebraic structure of the associated discrete (matrix) operators and their interrelations. We also stress the significant distinction between the discrete operators and the commonly adopted discretized versions of the continuous differential operations obtained via the method of finite differences in numerical analysis; the latter generally lack the desirable topological behaviour [42].

Figure 1.ȱȱAn example graph shown along with its edge–node incidence matrix, *Ȇ*, and cycle (edgeȬ face incidence) matrix, *ȍ* . Note that *M*=5, *Q*=4 and *R*=2. **Figure 1.** An example graph shown along with its edge–node incidence matrix, *Π*, and cycle (edge-face incidence) ma‐ trix, *Ω*. Note that *M*=5, *Q*=4 and *R*=2.

incidence matrix, in general, is not TU [30–34]. However, the edge–face incidence matrix is TU when each edge is included in exactly two faces that traverse the edge in opposite directions (e.g., a planar graph with a minimum cycle basis [30]). In this circumstance, total unimodularity of the edge–face incidence matrix stems from the fact that the face–edge incidence matrix is the edge–node incidence matrix of the dual graph [30]. Indeed, a TU constraint matrix (and integer constraints) guarantees that the solution of the related optimization problem (see also optimization problems (21) and (22) in the following) will be integer. Nonetheless, TU property has a further practical significance since the relaxed problem, obtained by neglecting the integer constraints, can also be solved using generic linear (not integer) programming solvers.ȱȱ As a final remark, some considerations on the *total unimodularity* (TU) property inherent to the matrix operators, which is extremely important in integer linear programming, are in order. We recall that a matrix is TU if the determinant of every submatrix is either zero or ±1. For any graph, the edge–node incidence matrix is TU. On the contrary, the face–edge incidence matrix, in general, is not TU [30–34]. However, the edge–face incidence matrix is TU when each edge

Ȭ

,{ , } *t*<sup>1</sup> *t<sup>Q</sup>*

Ȭ Ȭ

*Ȭ Ȭ*

2.2 Rigorous gradientȬbased formulation of the MChȬPhU problem

Ȭ

Ȭ

{ , , } *b*A<sup>1</sup> *b*A*<sup>Q</sup>*

is included in exactly two faces that traverse the edge in opposite directions (e.g., a planar graph with a minimum cycle basis [30]). In this circumstance, total unimodularity of the edge– face incidence matrix stems from the fact that the face–edge incidence matrix is the edge–node incidence matrix of the dual graph [30]. Indeed, a TU constraint matrix (and integer constraints) guarantees that the solution of the related optimization problem (see also optimization problems (21) and (22) in the following) will be integer. Nonetheless, TU property has a further practical significance since the relaxed problem, obtained by neglecting the integer constraints, can also be solved using generic linear (not integer) programming solvers.

## **2.2. Rigorous gradient-based formulation of the MCh-PhU problem**

Once the basic concepts of the discrete calculus and graph theory are presented, we are in the position to frame the formulation of the MCh-PhU problem on an appropriate mathematical playground. Let us consider a set of *Q single-look-complex* (SLC) SAR data acquired over a certain area of interest. One of them is assumed as the reference (master) image, with respect to which the images are properly coregistered. This set is characterized by the corresponding acquisition times {*t* 1 , …, *tQ*} and perpendicular baselines {*b*⊥1, …, *b*⊥*Q*}. Accordingly, for each coregistered SLC pair, a multilook *phase interferogram* (suitably depurated by the flat-Earth and topography contributions, by using a priori information about the topography and acquisition geometry) can be produced [37]; however, a common practice (within the multitemporal SB InSAR class) is first to identify a suitable small-baseline subset of the relevant multibaseline (temporal and spatial-perpendicular baselines) interferometric-pair set [6]. This is done to confine the effect of decorrelation phenomena associated with inherent angular and temporal electromagnetic backscattering variations [38]. Furthermore, a subset of common pixels of the *M* interferograms are then usually identified via the estimated *coherence* [37,38], so that only *P* pixels characterized by relatively high coherence values are singled out. In other words, the coherence index, providing a quantitative estimation of the decorrelation effects, permits to discriminate in favor of the "reliable" pixels.

The final aim is to reconstruct the absolute (i.e., not restricted in the principal [–π,π) interval) interferometric phases values from the wrapped (i.e., observed only in the principal [–π,π) interval) interferometric phase pertinent to *M* multichannel interferograms.

The problem we are interested in can be naturally framed on a discrete setting. Indeed, one possibility is to regard the discrete set of *P* selected (typically sparsely distributed) coherent pixels as a set of nodes, V <sup>B</sup>, in the Euclidean (*azimuth, range*) plane, and the set of *Q* SAR acquisitions representing a set of nodes, V <sup>A</sup>, in the Euclidean (*temporal-baseline, perpendicularbaseline*) plane [26]. Accordingly, a formal description of the problem at hand can be given in terms of a couple of abstract graphs: G<sup>A</sup> =(V <sup>A</sup>;ℰ<sup>A</sup> ) and G<sup>B</sup> =(V <sup>B</sup>;ℰ<sup>B</sup> ) where the corresponding edge sets ℰ<sup>A</sup> and ℰ<sup>B</sup> have to be properly defined. Accordingly, with *Q* and *M*, we denote the cardinality of V <sup>A</sup> and ℰ<sup>A</sup> , and with *P* and *N* the cardinality of V <sup>B</sup> and ℰ<sup>B</sup> , respectively. Note also that defining ℰ<sup>A</sup> pertains to the *M* interferometric pairs selection. Defining a meaningful edge set for a collection of nodes is now concerned since different criteria can be adopted to achieve it. The dimensionality of the ambient space in which the graph is embedded deserves some considerations. In this regard, we recall that a graph is called *planar* if it can be embedded in the plane [30–33]. Note also that a graph is not generally guaranteed to be planar, even if the

nodes are embedded in two dimensions. Since planarity is important for the workability of the implicated optimization procedure with powerful numerical solvers, a typical option to preserve graph planarity is resorting to the *Delaunay triangulation* in the Euclidean plane for establishing the edge set from nodes embedded in two dimensions [26]. Note that such an option specifically pertains to the solution strategy [26, 35]; nonetheless, our general formu‐ lation applies as well when different edge structures are adopted. Accordingly, once G<sup>A</sup> and G<sup>B</sup> have been somehow defined, the topological properties inherent to each graph are alge‐ braically captured by the related differential operators, which are summarized in Table 1.


**Table 1.** Adopted notation

First of all, we consider the absolute phase relevant to the multichannel SAR acquisition as a node variable pertinent to both the graphs GA and G<sup>B</sup> ; by using a matrix representation, this information can be conveniently arranged in a *Q* × *P* matrix *Φ* as follows:

$$\Phi = [\![\![\mathbf{q}^1, \dots, \mathbf{q}^P] \!] = \begin{bmatrix} \![\![\mathbf{q}\_1 \ ] \!] \\ \vdots \\ \![\![\mathbf{q}\_Q \ ] \!] \end{bmatrix} \tag{5}$$

where ∀ *p* ∈{ 1, 2, ..., *P*} **φ** *<sup>p</sup>* ∈ℝ*<sup>Q</sup>* encodes in a vectorial manner the *p*th node variable relevant to the graphs G<sup>A</sup> ; similarly, ∀*q* ∈{ 1, 2, ..., *Q*} **φ***<sup>q</sup>* ∈ℝ*<sup>P</sup>* encodes the *q*th node variable relevant to GB.

Widely adopted global gradient-based PhU approaches, which have historically been devel‐ oped for the single-channel case, generally consist in three processing steps [1, 29]. First, an estimation of the (wrapped) phase gradient is obtained; the estimated phase gradient is then suitably corrected (in terms of 2π multiples),and subsequently integrated to attain the unwrapped (absolute) phase.

Within the formulation of the MCh-PhU problem we concern [26], a twofold estimation of the discrete gradient field is carried out onto the considered two graphs GA and G<sup>B</sup>, as discussed in the following. The stack of the absolute interferometric phases relevant to the *M* (vectorized) interferograms can be formally represented through a *P* × *M* matrix denoted by

$$\begin{aligned} \boldsymbol{\Psi}' = [\boldsymbol{\updownarrow}\boldsymbol{\updownarrow}^1, \dots, \boldsymbol{\updownarrow}\boldsymbol{\updownarrow}^M] = \begin{bmatrix} \boldsymbol{\updownarrow}\boldsymbol{\updownarrow} \\ \vdots \\ \boldsymbol{\updownarrow}\boldsymbol{\updownarrow} \end{bmatrix} \end{aligned} \tag{6}$$

wherein the *P*-dimensional vector **ψ** *m* refers to the absolute phase field pertinent to the *m*th interferogram. Accordingly, *Ψ* is formally related to the absolute (unwrapped) phase matrix *Φ* via the discrete gradient operator *Π*<sup>A</sup> :

$$
\mathbf{M}^{\mathsf{T}} = \Pi\_{\mathsf{A}} \Phi \tag{7}
$$

Note also that

$$\mathbf{P}^{\mathsf{T}} = [\boldsymbol{\upuppi}^{\mathsf{T}}\_{1}, \dots, \boldsymbol{\upuppi}^{\mathsf{T}}\_{p}] = [\boldsymbol{\upuppi}\_{\mathsf{A}} \boldsymbol{\upup}^{\mathsf{I}}, \dots, \boldsymbol{\upuppi}\_{\mathsf{A}} \boldsymbol{\upup}^{p}] \tag{8}$$

By applying the discrete gradient operator *Π*<sup>B</sup> to the absolute phase of each interferometric pair, we obtain

1 B B B [ , , ] *<sup>M</sup> Π Π Π* = **ψ ψ** K (9)

As a result, by using Eq. (7), we get

$$
\Pi\_{\mathsf{B}} \mathsf{\varPi}\_{\mathsf{B}} = \Pi\_{\mathsf{B}} \mathsf{\varPhi}^{\mathsf{T}} \Pi\_{\mathsf{A}}^{\mathsf{T}} \tag{10}
$$

Second, we consider the (wrapped) interferometric phase that is uniquely defined only in the principal value range since it is obtained as the phase of a complex function. Hence, it is convenient to formally introduce the non-injective (modulo-2π) *wrapping operator W : φ ∈ ℝ → mod (φ + π, 2π) − π ∈ − π, π).* It should be noted that the following trivial identi‐ ties hold:

$$\begin{aligned} \mathcal{W}(\mathcal{W}(A) \pm \mathcal{W}(\mathcal{B})) &= \mathcal{W}(A \pm \mathcal{B}) = \mathcal{W}(\mathcal{W}(A) \pm \mathcal{B}) = \mathcal{W}(A \pm \mathcal{W}(\mathcal{B})) & \mathbf{a} \\ \mathcal{W}(A) &= A + 2\pi \mathbf{Z} & \mathbf{b} \\ \hline \end{aligned} \tag{11}$$

where *A* and *B* represent two generic matrices and *Z* is a suitable integer matrix. Given the stack of the unknown absolute (unwrapped) interferometric-phases *Ψ*, the corresponding stack of the wrapped phases *Ψ***˜** can be conveniently expressed in terms of the *wrapped discrete gradient of* (*wrapped*) *observed phase* as follows:

$$\tilde{\mathbf{P}}^{\mathsf{T}} = [\mathcal{W}(\varPi\_{\mathcal{A}} \mathcal{W}(\mathsf{q}^{1})), \dots, \mathcal{W}(\varPi\_{\mathcal{A}} \mathcal{W}(\mathsf{q}^{p}))] = \mathcal{W}(\varPi\_{\mathcal{A}} \mathcal{W}(\mathsf{q})) = \mathcal{W}(\varPi\_{\mathcal{A}} \boldsymbol{\Phi}) \tag{12}$$

where we have exploited Eq. (11b). Note also that the *p*th column of *Ψ***˜** Τ , i.e., **ψ˜** *p* T =*W* (*Π*A*W* (**φ** *p* )), can be regarded as an estimate of the absolute-phase discrete gradient on the graph G<sup>A</sup>. It is worth remarking that the observed (multilook) interferometric phase can however be corrupted by noise [39–41], which is taken into account by considering an additive phase noise term *D*. Accordingly, by using Eq. (11a), we get

$$\tilde{\mathbf{P}}^{\mathsf{T}} = \mathsf{W}(\mathsf{W}(\varPi\_{\mathsf{A}}\varPhi) + \mathsf{D}) = \mathsf{W}(\varPi\_{\mathsf{A}}\varPhi + \mathsf{D})\tag{13}$$

More specifically, whenever a possible spatial filtering (e.g., conventional multilooking followed by a noise-filtering step [42]) is *independently* applied to each SAR interferometric data pair, the resulting term *D* in Eq. (13) implies that the phase interferograms **ψ˜** *m* , with *m*∈{ 1, 2, ..., *M* }, are no more fully time consistent (in the sense of [26, 43, 44]). To clarify this point, we observe that by using Eq. (13), it turns out that

$$\mathcal{W}\left(\mathbf{D}\_{A}^{T}\mathbf{P}^{T}\right) = \mathcal{W}\left(\mathbf{D}\_{A}^{T}\left(\mathbf{D}\_{A} + \mathbf{D} + 2\pi\mathbf{Z}\right)\right) = \mathcal{W}\left(\mathbf{D}\_{A}^{T}\mathbf{T}\_{A}\mathbf{D} + \mathbf{D}\_{A}^{T}\mathbf{D}\right) = \mathcal{W}\left(\mathbf{D}\_{A}^{T}\mathbf{D}\right) \neq \mathbf{0} \tag{14}$$

where we have used Eq. (11b) with *A*=*Π*A*Φ* + *D* and noted that *Ω*<sup>A</sup> <sup>T</sup> *Z* is an integer matrix and *Ω*<sup>A</sup> <sup>T</sup>*Π*<sup>A</sup> =**0** (according to Eq. (3)). Eq. (14) reads as "the wrapped discrete curl of the interfero‐ metric phase on G<sup>A</sup> (i.e., pertinent to the 'temporal' domain) is generally different from zero"; it formally expresses the (modulo-2π) *cyclic inconsistency* of the multichannel interferometric phase inherent to independently filtered SAR interferograms. Note also that Eq. (14) repre‐ sents, within our framework, the generalization (to a wider class of discrete settings) of the "phase triangularity" condition in ref. [44], capturing the underlying structure of the problem within a suitable matrix formalism. With reference to the *m*th interferometric pair, the

estimated absolute interferometric-phase gradient on the graph G<sup>B</sup> is then obtained by *wrapping the discrete gradient of* (*wrapped*) *interferometric-phase field*: **g** *m* <sup>=</sup>*<sup>W</sup>* (*Π*B**ψ˜** *m* ). Thus, by stacking the so-obtained absolute phase gradient estimations, we get the *N* × *M* matrix *G* = **g** 1 , **g** 2 , …, **g** *<sup>M</sup>* , where

$$\begin{array}{|c|c|c|}\hline \cr & \cr & \mathbf{G} = \mathbf{W(H\_{z}\tilde{\mathbf{W}})} \\ \hline \cr & \cr & \cr & \mathbf{G} \\ \hline \cr & \cr & \mathbf{G} & \mathbf{G} \\ \hline \end{array} \qquad \begin{array}{|c|c|}\hline \mathbf{G=W(H\_{z}\tilde{\mathbf{W}})} \\ \hline \cr & \cr & \mathbf{G} \\ \hline \cr & \mathbf{G} & \mathbf{G} \\ \hline \cr & \mathbf{G} & \mathbf{G} \\ \hline \end{array} \qquad \begin{array}{|c|c|}\hline \cr & \cr & \mathbf{G} \\ \hline \cr & \mathbf{G} & \mathbf{G} \\ \hline \cr & \mathbf{G} & \mathbf{G} \\ \hline \end{array} \qquad \begin{array}{|c|c|}\hline \cr & \mathbf{G} \\ \hline \cr & \mathbf{G} \\ \hline \cr & \mathbf{G} \\ \hline \end{array} \qquad \begin{array}{|c|c|}\hline \cr & \mathbf{G} \\ \hline \cr & \mathbf{G} \\ \hline \cr & \mathbf{G} \\ \hline \end{array}$$

$$\mathbf{G} = \mathcal{W} \left( \mathbf{Z} \mathbf{I}\_{\mathbf{B}} \tilde{\mathbf{P}}^{\top} \right) = \mathcal{W} \left( \mathbf{Z} \mathbf{I}\_{\mathbf{B}} \mathcal{W} \left( \left[ \mathbf{Z} \mathbf{I}\_{\mathbf{A}} \boldsymbol{\Phi} \right]^{\top} + \mathbf{D}^{\top} \right) \right) \tag{16}$$

From Eq. (16), by using Eq. (11b), we get

$$\mathbf{G} = \mathcal{W} \left( \mathbf{T}\_{\mathrm{B}} \boldsymbol{\Phi}^{\mathrm{T}} \mathbf{T}\_{\mathrm{A}}^{\mathrm{T}} + \mathbf{T}\_{\mathrm{B}} \mathbf{D}^{\mathrm{T}} + 2\pi \,\mathbf{T}\_{\mathrm{B}} \mathbf{Z} \right) = \mathcal{W} \left( \mathbf{T}\_{\mathrm{B}} \boldsymbol{\Phi}^{\mathrm{T}} \mathbf{T}\_{\mathrm{A}}^{\mathrm{T}} + \mathbf{T}\_{\mathrm{B}} \mathbf{D}^{\mathrm{T}} \right) \tag{17}$$

where in the last equality we have noted that *Π*B*Z* is also an integer matrix. It should be emphasized that, under the assumption *D* =*o*, the equality between Eqs. (10) and (17) holds only up to an integer matrix multiplied by 2*π*.

#### **2.3. The MCh-PhU problem as constrained optimization**

In this Section, the nonlinear inversion MCh-PhU problem is reformulated as a (nonlinear) constrained optimization problem. According to the presented general formulation, we introduce in the following the MCh-PhU problem as the solution of the following matrix equation:

$$\left\lfloor \begin{array}{c} \square \\ \square \end{array} \right\rfloor \square \left\lfloor \begin{array}{c} \square \\ \square \end{array} \right\rfloor \square \left\langle \begin{array}{c} \square \\ \square \end{array} \right\rfloor \left\langle \square \\ \square \end{array} \right\rfloor \left\langle \begin{array}{c} \square \\ \square \end{array} \right\rangle \left\langle \square \\ \square \end{array} \right\rangle \left\langle \square \\ \square \left\{ \begin{array}{c} \square \\ \square \end{array} \right\rangle \left\langle \square \\ \square \end{array} \right\rangle \left\langle \square \\ \square \square \right\rangle \left\langle \square \\ \square \square \end{array} \right\rangle \left\langle \square \begin{array}{c} \square \\ \square \square \\ \square \square \end{array} \right\rangle \left\langle \square \\ \square \square \square \square \end{array} \right\rangle \begin{array}{c} \square \\ \square \square \\ \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square \square$$

where the columns of *G* represent the interferometric-phase *pseudo*-gradients estimated from the observed phase, and *K* is an (unknown) *N* × *M* integer matrix, whose columns represent the corresponding (2π-normalized) corrections to be added to the (wrapped) interferometricphase *pseudo*-gradients in order to recover the absolute interferometric-phase discrete gradi‐ ents. It is worth noting that the term pseudo-gradient is used here to emphasize that the integration of the estimated gradient is path-dependent (non-conservative behavior); the term *residues* [1] is also typically used to connote the inconsistency of the estimated phase gradient. As a matter of fact, matrix equation (Eq. (18)) describes an *ill-posed* problem, in which the data *G* generally do not constrain sufficiently the problem to get a unique solution. Additional suitable constraints and *a priori* assumption have, thus, to be introduced to solve the problem. First, for restoring the *cyclic consistency* (see Section 2.1) of the estimated pseudo-gradients pertinent to the graphs GB and G<sup>A</sup>, two corresponding sets of (equality) constraints have to be enforced, respectively. More specifically, pre-multiplying both sides of Eq. (18) by *Ω*<sup>B</sup> T and taking into account Eq. (3), we obtain

( ) <sup>T</sup> <sup>B</sup> *Ω G K 0* - = 2p(19)

Similarly, by premultiplying both sides of the transposed version of Eq. (18) by *Ω*<sup>A</sup> T and taking also into account Eq. (3), we obtain

$$\mathbf{Q}\_A^T \left( \mathbf{G} - \boldsymbol{\Pi}\_\mathbf{B} \mathbf{D}^T - 2\boldsymbol{\pi} \, \mathbf{K} \right)^\mathrm{T} = \boldsymbol{\theta} \tag{20}$$

Constraints stated by Eq. (19) imply that the columns of *G* −2*π K* must lie in the null-space of *Ω*<sup>B</sup> T . Since the matrix *Π*<sup>B</sup> represents a basis to span the null-space of *Ω*<sup>B</sup> T (see Eq. (3)), we may then write *G* −2*π K* = *Π*B*X* , where *X* is a new variable. Accordingly, the corrected pseudogradients stack *G* −2*π K* is enforced to be a stack of discrete gradients, which can thus be unambiguously integrated. Similarly, Eq. (20) implies *G* −*Π*B*D* T −2*π K* T =*Π*A**Y**. As a result, the two sets of constraints, stated by Eqs. (19) and (20), guarantee that the solution of the problem is effective in preserving the *cyclic consistency* (curl-free) property of the corrected gradients pertaining to the graphs GB and G<sup>A</sup>, respectively. As a matter of fact, the solution of Eq. (18) cannot be determined by using the two sets of constraints (Eqs. (19) and (20)) only; thus, the inverse problem must be first regularized [45]. The minimum-norm methods search for a global solution that minimizes a generalized error-norm associated with an optimality criterion, so incorporating prior information about the behavior of the solution [1]. Accord‐ ingly, we resort to a regularization approach using *l*<sup>1</sup> -norm minimization in weighted version, as a specific case of *l*<sup>p</sup> -norm general formulation. Formally, the MCh-PhU problem may be then formulated as a constrained optimization problem for the field of integer corrections:

$$\hat{\mathbf{K}} = \operatorname\*{arg\,min}\_{\mathbf{K} \in \mathbb{Z}^{N \times M}} \left\| \mathbf{K} \right\|\_{1, \mathcal{C}} \tag{21}$$

subject to

$$\begin{cases} \boldsymbol{\mathfrak{Q}}\_{\text{A}}^{\mathrm{T}} \boldsymbol{K}^{\mathrm{T}} = \boldsymbol{\mathfrak{Q}}\_{\text{A}}^{\mathrm{T}} \left[ \boldsymbol{G} - \boldsymbol{\varPi}\_{\text{B}} \boldsymbol{\mathsf{D}}^{\mathrm{T}} \right]^{\mathrm{T}} \left( 2\boldsymbol{\pi} \right)^{-1} \\\ \boldsymbol{\mathfrak{Q}}\_{\text{B}}^{\mathrm{T}} \boldsymbol{K} = \boldsymbol{\mathfrak{Q}}\_{\text{B}}^{\mathrm{T}} \boldsymbol{G} \left( 2\boldsymbol{\pi} \right)^{-1} \end{cases} \tag{22}$$

wherein

$$\left\|\mathbf{K}\right\|\_{1,\mathcal{C}} = \sum\_{m=1}^{M} \sum\_{n=1}^{N} c\_{nm} \left| k\_{nm} \right| \tag{23}$$

represents the weighted *l*<sup>1</sup> -norm [46] of the matrix *K*, *C* = *cnm <sup>N</sup>* <sup>×</sup>*<sup>M</sup>* denotes a suitable weighting matrix, and ℤ indicates the field of integer numbers. As far as the existence of an integer solution for Eqs. (21) and (22) is concerned, it should be noted that the considerations at the end of Section 2.1 apply. Since the first matrix equation in Eq. (22) includes a generally not null (unwanted) term *Ω*<sup>A</sup> <sup>T</sup>*D* , its fulfillment deserves further discussion. Although the evaluation of *W* (*Ω*<sup>A</sup> <sup>T</sup>*D*) can be obtained according to Eq. (14), however, a full estimation for *Ω*<sup>A</sup> <sup>T</sup>*D* is generally not a simple task. Further discussion is provided in Section 3. The solution of the optimization problem (Eqs. (21) and (22)) is also referred to as the minimum weighted discontinuity solution (in a weighted *l*<sup>1</sup> -norm sense) [1, 23]. As a matter of fact, finding the global minimum point of the problem stated by Eqs. (21) and (22) for an arbitrary pair of graphs is, in general, a difficult task. A suboptimal strategy aimed at solving Eqs. (21) and (22) consists in adopting a two-stage approach. This is, in particular, the solution strategy implemented through the *extended minimum cost flow* (EMCF) technique [26], in which the edge structure of each considered graph is usually defined via a *Delaunay* triangulation in the Euclidean plane, to take advantage from efficient solvers for minimum cost flow (MCF) problems [47, 49, 50]. We remark that the distinctive characteristic of the EMCF approach is the extensive use of the computationally efficient MCF method. Moreover, a dual-level parallel model for EMCF has also been proposed in refs. [35] and [36]. Moreover, different approaches toward full 3D phase unwrapping have recently been proposed in refs. [63] and [64].

## **3. Noise-filtering of multichannel SAR interferograms**

In this Section, we review the basic concepts concerning the filtering of noise that corrupts a stack of multitemporal SAR interferograms. First, the noise-filtering operation for singlechannel multilook interferograms is discussed; subsequently, the general framework of the multichannel noise-filtering (MC-NF) approach, which is intimately connected with the problem of multichannel phase unwrapping, is described.

#### **3.1. Decorrelation noise in SAR interferograms**

In order to introduce the problem at hand, let us first consider one single-channel SAR interferogram obtained starting from two SAR (synthetic aperture radar) images, namely, *i* 1 and *i* 2 , acquired (over the same scene on Earth) at two different times, namely, *t* 1 and *t* 2 , respectively. The two SAR images can be represented via two complex-valued signals, say *i*(*x*, *r*) and *i* 2 (*x*, *r*), with *x* and *r* denoting the two independent spatial variables (with respect

to azimuth and range direction, respectively) in the radar geometry. The two complex signals can be expressed as follows [29, 51]:

$$\begin{array}{c} \mathrm{i}\_{1}(\mathbf{x},r) = \mathbf{y}\_{1}(\mathbf{x},r)e^{-j\frac{4\pi}{\lambda}r} + \mathrm{n}\_{1}(\mathbf{x},r) \\\\ \mathrm{i}\_{2}(\mathbf{x},r) = \sum\_{\begin{subarray}{c} \mathbf{i}\_{2}(\mathbf{x},r) = \mathbf{y}\_{2}(\mathbf{x},r) \\ \mathbf{x} \leq \mathbf{y}\_{1} \leq \mathbf{y}\_{2} \end{subarray}} \mathrm{i}\_{2}(\mathbf{x},r) = \mathrm{y}\_{2}(\mathbf{x},r)e^{-j\frac{4\pi}{\lambda}(r+\delta r)} + \mathrm{n}\_{2}(\mathbf{x},r) \\\\ \mathrm{i}\_{3}(\mathbf{x},r) = \sum\_{\begin{subarray}{c} \mathbf{i}\_{2}(\mathbf{x},r) = \mathbf{y}\_{1} \leq \mathbf{y}\_{2} \end{subarray}} \left( \begin{array}{c} \mathbf{x},r \\ \mathbf{x} \end{array} \right) \left( \begin{array}{c} \mathbf{x},r \\ \mathbf{x} \end{array} \right) \left( \begin{array}{c} \mathbf{x},r \\ \mathbf{x} \end{array} \right) \end{array} \tag{24}$$

where *δr* is the sensor-to-target slant range difference at time *t* 2 with respect to time *t* 1 , and *γ*1 (*x*, *r*) and *γ*<sup>2</sup> (*x*, *r*) are the corresponding (complex-valued) reflectivity functions of the illuminated scene at time *t* 1 and *t* 2 , respectively. Furthermore, the two additive (noise) contributions *n*<sup>1</sup> (*x*, *r*) and *n*<sup>2</sup> (*x*, *r*) describe random quantities that are included in Eq. (24). As a result of these noise terms and of the intrinsic random nature of the two images reflectivity functions, when the two SAR images are interfered to form a so-called interferogram, i.e., when their phase difference *ψ* is extracted, the interferometric phase will be noisy. An important parameter influencing the quality of the retrieved interferometric phase is the (complex) *crosscorrelation* factor between the two involved SAR images, which is typically defined as

$$\mathcal{X} = \frac{E[\dot{\boldsymbol{i}}\_1 \cdot \dot{\boldsymbol{i}}\_2^\*]}{\sqrt{E\left[\left|\dot{\boldsymbol{i}}\_1\right|^2\right] \cdot E\left[\left|\dot{\boldsymbol{i}}\_2\right|^2\right]}} = \rho \,\mathrm{e}^{j\nu} \tag{25}$$

where *ρ* ∈ 0, 1 , *ψ* ∈ −*π*, *π*), and the asterisk denotes the conjugate complex value. Notewor‐ thy, the cross-correlation factor (Eq. (25)) is a complex-valued term that can be decomposed in terms of amplitude *ρ* (i.e., *ρ* = |*χ* |) and phase *ψ*. For interferometric SAR images, *χ* can be evaluated by performing spatial averaging (known as multilooking) operations on a *statisti‐ cally homogeneous* area. Indeed, the symbol *E* in Eq. (25), which is representative of the statistical expectation operation [52], can then be replaced by the spatial averaging operation. The amplitude factor *ρ*, which is known to as *coherence*, accounts for the similarity between the two SAR images, whereas *ψ* is the multilook interferometric phase. A value for the coherence that approaches zero is representative of an uncorrelated scene, whereas coherence value that is close to unity corresponds to a noise-free interferogram.

There are several causes that are responsible for coherence decrease. As matter of fact, the cross-correlation factor in Eq. (25) depends on different noise sources, and it can be conven‐ iently factorized as follows [51]:

$$\mathcal{X} = \mathcal{X}\_{\text{the}} \cdot \mathcal{X}\_{\text{term}} \cdot \mathcal{X}\_{\text{spa}} \cdot \mathcal{X}\_{\text{dop}} \cdot \mathcal{X}\_{\text{mis}} \cdot \mathcal{X}\_{\text{vol}} \cdot e^{\mathcal{W}} \tag{26}$$

#### where


The multilook operation, leading to the multilook phase *ψ* in Eq. (25), reduces the level of noise corrupting interferograms, although this is paid in terms of a reduction of spatial resolution of interferograms. Multilook interferometric phase can be described by using a random quantity and, accordingly, it can be characterized via the knowledge of its probability density function. It has been shown in literature [40, 41, 54] that the *probability density function* (pdf) of an *L*-multilook interferometric phase (with *L* being the number of averaging samples in the averaging window used for the estimation of the statistical average operation involved in the calculation of Eq. (25)) can be given in terms of a Gauss *hypergeometric* function:

$$p\_L(\boldsymbol{\nu}) = \frac{\Gamma\left(L + \frac{1}{2}\right) (1 - \rho^2)^L \beta}{2\sqrt{\pi} \,\Gamma\left(L\right) \left(1 - \beta^2\right)^{L + \frac{1}{2}}} + \frac{(1 - \rho^2)^L}{2\pi} \,\_2F\_1\left[L, 1; \frac{1}{2}; \beta^2\right] \tag{27}$$

where *ψ* ∈ −*π*, *π*), *β* =*ρcos*(*ψ*−*ψ*<sup>0</sup> ), *ρ* represents the coherence, *F*<sup>2</sup> <sup>1</sup> denotes the Gauss hyper‐ geometric function, Γ(⋅) is the *gamma* function, and *ψ*<sup>0</sup> represents the expected "true" value of the interferometric phase. The peak of the distribution is located at *ψ*=*ψ*<sup>0</sup> .

The pdf in Eq. (27) is sketched in Figure 2a for different values of *L* and in Figure 2b for different values of the *ρ*, with *ψ*<sup>0</sup> =0. By observing Figure 2a, it is clear that pdfs become narrower as the number of looks *L* increases (as expected). This finding is extremely important because it demonstrates that the interferometric phase may be thought to be corrupted by an additive noise random signal, namely, *v*, that has the same pdf as in Eq. (27) but with a zero-mean expected value, i.e., we may assume as valid the following additive model for the interfero‐ metric noise [54]: *ψ*=*ψ*<sup>0</sup> + *v*. To further investigate about the statistics of multilook interfero‐ grams, we can observe that the validity of Eq. (27) is only restricted to the [–π,π) interval. However, this restriction does not apply when the phase signal is directly derived in the complex plane instead of the real plane. In the works of Lopez (2003) [55] and Lopez and Pottier (2007) [40], a comprehensive analytical derivation of the noise statistics in the complex plane is derived. Nonetheless, the *Cramér–Rao* bound for the standard deviation of multilook phase is given by [59]

$$
\sigma\_v = \frac{1}{\sqrt{2L}} \frac{\sqrt{1 - \rho^2}}{\rho} \tag{28}
$$

that shows that standard deviation depends on the coherence *ρ* and multilook factor *L.* Note that the phase standard deviation approaches the limit (Eq. (28)) asymptotically as the number of looks increases.

Figure 2. Probability density function of the interferometric phase [rad]: (a) for different values of the number of looks *L* (1—black line, 2—red line, 5—blue line, 10—green line, 20—orange line), with ; (b) for different values of the correlation coefficient (0.1—black line, 0.2—red line, 0.5 blue line, 0.7—green line, 0.8—orange line) with *L* = 4. **Figure 2.** Probability density function of the interferometric phase *ψ* [rad]: (a) for different values of the number of looks *L* (1—black line, 2—red line, 5—blue line, 10—green line, 20—orange line), with *ρ* =0.7 ; (b) for different values of the correlation coefficient *ρ* (0.1—black line, 0.2—red line, 0.5—blue line, 0.7—green line, 0.8—orange line) with *L* = 4.

$$
\sigma\_{\rm v} = \frac{1}{\sqrt{2L}} \frac{\sqrt{1 - \rho^2}}{\rho} \tag{28}
$$

that shows that standard deviation depends on the coherence and multilook factor *L*.

‐ ‐

‐ ‐

‐ ‐

### **3.2. Single-channel noise-filtering approaches for multilook interferograms**

In order to mitigate the effects of decorrelation noise artifacts affecting SAR interferograms, several noise-filtering techniques, mostly working on single-channel data, have been proposed in literature over the years [42, 54, 56, 57]. As shown in previous Section 3.1, the statistics of multilook interferograms can be characterized via a probability density function expressible in closed form (Eq. (27)), and the noise standard deviation generally depends upon the *coherence ρ* and the *number of looks L* [see also Eq. (28)]. Three different multilook interferograms, which are characterized by the same perpendicular baseline (of about 100 m), have been obtained by using three SAR sensors working at the different (C, X, and L) bands of the microwave region and are depicted in Figure 3. As it is evident from Figure 3, L-band inter‐ ferograms are less affected by noise than the ones pertinent to C and/or X-band.

**Figure 3.** Multilook interferogram computed by using different SAR data pairs: (a) July 11, 2011–August 16, 2011, Xband Cosmo-SkyMed (CSK); (b) September 15, 2004–October 29, 2004, C-band ASAR/ENVISAT; (c) July 30, 2007–Sep‐ tember 14, 2007, L-band ALOS/PALSAR-1.

It should be emphasized that coherence and noise levels can also significantly change from one SAR system to another depending on the operational wavelength.Multilook processing has been proved to be effective for noise reduction, but this is paid in terms of a decrease of the image spatial resolution. Noise filtering constitutes a preliminary step before phase unwrapping. Indeed, the multilook operation usually involves an averaging on neighboring SAR pixels, hence reducing the spatial resolution of the interferograms. Several algorithms have been proposed in literature. The most commonly used noise filter is the *boxcar* filter applied in the complex plane. Another frequently used option is provided by the *Golstein*'s frequency-domain algorithm [42], which is an empirical approach proposed for geophysical applications. In this case, a complex interferogram (amplitude and phase) is segmented into overlapping rectangular patches and for each patch the relevant power spectrum *Z* is com‐ puted.

**Figure 4.** Multilook interferogram relevant to the SAR data pair July 11, 2011–August 16, 2011, X-band Cosmo-SkyMed (CSK): (a) Original, (b–d) Goldstein filtering, with (b) *α* =0.25, (c) *α* =0.5, and (d) *α* =1.0.

The response of the *Golstein*'s filter is then computed from the power spectrum as follows:

$$\mathcal{H}(\xi,\eta) = \left| Z(\xi,\eta) \right|^a \tag{29}$$

where *ξ* and *η* denote the relevant spectral variables, respectively. Notice that the filtering effect vanishes when *α* =0 ; conversely, the filtering effect is more pronounced as the parameter *α* approaches unity. We show in Figure 4 the result of the application of the Goldstein's filter to a multilook interferogram, relevant to the Mt. Etna (Italy) volcano, obtained by using the *Cosmo-Skymed* sensor of the Italian Space Agency (ASI). Specifically, different values of the filtering parameter *α* have been considered in Figure 4. The limited effectiveness of the filtering capabilities of Goldstein approach is evident from the result depicted in Figure 4. A modifi‐ cation of the Goldstein filter that relies on an adaptive choice of the filtering factor *α* (which depends on the spatial coherence *ρ*) has also (more recently) been proposed by Baran in 2003 [58]. Other filters, such as the *median* filter [59] and the two-dimensional *Gaussian* filter, are also used to reduce noise while performing multilooking operations. It is worth noting that boxcar and Goldstein filters do not adapt to the direction of the fringes because these filters are operated in a square window. In order to overcome such a limitation, Lee et al. 1998 [54] then proposed a self-adaptive filter based on local gradient slope estimation that is able to improve noise-filtering performance by exploiting directional characteristics of an InSAR interferogram. Several adaptations and relevant improvements of the *Lee filter* have subse‐ quently proposed in literature over the recent years [56, 57], most of them based on the exploitation of the intrinsic directional behavior of InSAR interferograms. In fact, compared with the fringe phase and gradient, the fringe direction variation is gently, thus making it possible to use fringe direction to guide interferogram filtering.

#### **3.3. The multichannel noise-filtering (MCh-NF) algorithm**

The noise-filtering methods discussed in the previous Section have historically been developed to analyse and filter out the noise affecting single interferograms, separately, thus without taking into account their mutual temporal relationships. A multichannel noise-filtering problem arises when a stack of SAR interferograms need to be jointly exploited. In this case, it is profitable to develop/use noise-filtering approaches that not only exploit spatial/frequency information but can also take into account temporal relationships among available multichan‐ nel interferograms, in order to identify and filter out the noise affecting the whole interfero‐ metric data stack in the more reliable way as possible. A specific multichannel noise-filtering (MCh-NF) method [43], which is based on using a stack of time-redundant multilook inter‐ ferograms, is described in this Section. The MCh-NF method is here described by adopting the same rigorous formalism and terminology used for the topological description of multichannel phase unwrapping problem presented in Section 2. According to the adopted symbolism, let us consider *Q* SAR images and let *M* be the number of multilook interferograms characterized by small perpendicular and temporal baselines.

The resulting interferometric data stack of the *M* (wrapped) small-baseline multilook inter‐ ferograms can be expressed as *Ψ***˜** Τ <sup>=</sup> **ψ˜** 1 T , …, **ψ˜** *P* T ; thus, the *M*-dimensional vector **ψ˜** *p* T described the (vectorized) multichannel interferometic-phase pertinent to the *p*th pixel, with *p* ∈{1, …, *P*} and *P* denoting the number of coherent pixels common to all interferograms. In particular, *Ψ***˜** can be expressed in terms of discrete gradient *Π*<sup>A</sup> , according to Eq. (13), as:

$$\begin{array}{c} \begin{array}{c} \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \\ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \end{array} \circ \begin{array}{c} \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \\ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{array} \\ \end{array} \circ \begin{array}{c} \begin{array}{c} \\ \end{$$

wherein *Φ* represents the (unknown) phases associated with the available SAR images, and the matrix *D* describes the additive noise-term that corrupts the stack of interferograms. The noise term should be estimated and properly filter out from the generated interferograms. As discussed in Section 3.2, the term *D* arises since both a multilook operation and a noise-filtering procedure are typically applied to each single interferogram, separately. Both these operations are independently carried out on each single interferogram; hence, they are not necessarily time consistent. The fact that the interferograms are not fully time consistent can be formally expressed, according to Eq. (14), in terms of discrete curl *Ω*<sup>A</sup> T , in the form:

$$\mathcal{W}\left(\mathbf{\mathcal{Q}}\_A^T \tilde{\mathbf{P}}^T\right) = \mathcal{W}\left(\mathbf{\mathcal{Q}}\_A^T \mathbf{D}\right) \neq \mathbf{0} \tag{31}$$

which represents the *topological* generalization of the phase-triangularity condition exploited by the *SqueeSAR* technique [44]. Therefore, the *multichannel noise-filtering* (MCh-NF) approach suitably addresses the temporal inconsistencies inherent in the time-redundant multilook interferograms, which can be mathematically described in terms of the (modulo-2π) *cyclic inconsistency* of the multichannel interferometric phases [see Eq. (31)]. More specifically, MCh-NF is based on the solution of the a *nonlinear* optimization problem, as detailed in the following. First, ∀ *p* ∈{1, …, *P*}, the *Q*-dimensional vector (*Q* is the number of SAR acquisitions) repre‐ senting the (unknown) wrapped phases **Φ˜** *p* =*W* (**Φ***<sup>p</sup>* ) is estimated as follows:

$$\hat{\tilde{\mathbf{D}}}^{\boldsymbol{\nu}} = \underset{\tilde{\Phi}^{\boldsymbol{\nu}} \in \mathbb{R}^{Q}}{\arg\max} \quad \left| \overline{\mathcal{L}}\_{\boldsymbol{\nu}} \circ e^{\boldsymbol{\beta}(\tilde{\Psi}\_{\boldsymbol{\nu}}^{\boldsymbol{\Gamma}} - \mathcal{W}(\boldsymbol{H}\_{\boldsymbol{\Lambda}}\bar{\Phi}^{\boldsymbol{\nu}}))} \right| \tag{32}$$

where *j* = −1 denotes the imaginary unit, *P* denotes the number of coherent pixels to which the noise-filtering procedure is applied, the symbol represents the *Hadamard* product, and **ζp** ͟ = *ζ* ͟ *p* 1 , ⋯, *ζ* ͟ *p M T* is an *M*-dimensional *normalized* weighting vector representing our confidence on the quality of the exploited *M* (small-baseline) interferometric phases pertinent to the *p*th pixel, with

$$\overline{\mathcal{L}}\_p^m = \frac{\mathcal{L}\_p^m}{\sum\_{h=1}^M \mathcal{L}\_p^h} \tag{33}$$

wherein the generic elements *ζ<sup>p</sup> m* can be related to the *spatial coherence* as detailed after. Subsequently, these estimated vectors **Φ˜ ^** *p* are used to reconstruct a new (noise-filtered) stack of time-consistent interferograms *Ψ***˜ ^** Τ <sup>=</sup> **ψ˜ ^** 1 T , …, **ψ˜ ^** *P* T , where **Ψ˜ ^** *p T* <sup>=</sup>*<sup>W</sup>* (*Π*A**Φ˜ ^** *p* ) and *p* ∈{1, …, *P*}. We emphasize that, according to Eq. (32), the MCh-NF technique is based on searching for the (unknown) wrapped-phase vector **Φ˜** *<sup>p</sup>* ∈ℝ*<sup>Q</sup>* that minimizes the (weighted) circular variance of the random (phase) vector representative of the phase difference, **Ψ˜** *p T* <sup>−</sup>*<sup>W</sup>* (*Π*A**Φ˜** *p* ), between the "original" and the "reconstructed" interferograms.

The evaluation of the weights for the optimization problem in Eq. (32) is now addressed. Let *Θ* ͜ = Θ ͜ *i*, *j* be a matrix description for a generic 2-D phase map, whose corresponding vectorized representation is provided by the *P*-dimensional vector **Θ**. Each pixel of the phase map is identified by discrete range and azimuth coordinates, denoted by *i* and *j,* respectively. Accordingly, each pair (*i*, *j*) is uniquely associated with a monodimension index *p* ∈{1, …, *P*} identifying an element of the vector *Θ*. The *spatial coherence* relevant to **Θ** (i.e., Θ ͜ ) evaluated around the pixel (*i*, *j*) (associated with the index *p*) is defined as

$$\mathcal{L}\_p\left(\Theta\right) = \frac{1}{\left(2L\_R + 1\right)\left(2L\_A + 1\right)} \left| \sum\_{l=-L\_R}^{L\_R} \sum\_{k=-L\_A}^{L\_A} \exp\left[j\left|\bar{\Theta}\_{i+l, j+k}\right|\right] \right| \tag{34}$$

where 2*L<sup>R</sup>* + 1 and 2*L<sup>A</sup>* + 1 are the number of azimuth and range pixels within the used boxcar averaging window, which is centred around the generic pixel identified by the discrete range and azimuth coordinates, *i* and *j*, respectively. In particular, the *m*th weight *ζ<sup>p</sup> m* is expressed, according to Eq. (34), in terms of *spatial coherence* relevant to the (vectorized) interferograms **ψ˜** *m* and evaluated around the pixel associated with the index *p*, in the functional form *ζp m* =*ζ<sup>p</sup>* (**ψ˜** *m* ), ∀*m*∈{1, …, *M* }. Therefore, **ζ***<sup>p</sup>* = *ζ<sup>p</sup>* 1 , ⋯, *ζ<sup>p</sup> M* T can be evaluated in terms of the spatial coherence directly from the stack of *M* multilook interferograms *Ψ***˜** <sup>=</sup> **ψ˜** 1 , …, **ψ˜** *<sup>M</sup>* . As experimentally demonstrated in ref. [43], the "reconstructed" interferograms with MCh-NF are significantly less affected by noise than the original ones. However, a group of the reconstructed interferograms, although limited, can exhibit spatial coherence values smaller than the ones relevant to the corresponding original interferograms, thus implying that a partial corruption of the spatial coherence occurs during the minimization procedure. In particular, it happens in correspondence to interferograms that were originally significantly coherent, and this is due to the fact that the strategy in Eq. (32) tends to "inject" coherence from very coherent to incoherent interferograms, by exploiting the time redundancy of the small

baseline data pairs. Accordingly, in order to also preserve the spatial coherence of the very coherent interferograms, a simple nonlinear combination between the original and the reconstructed interferograms is carried out, thus further increasing the phase quality of the whole set of *M* reconstructed interferograms. In particular, the two sets of interferograms are combined through the following (wrapped) weighted averaging operation:

$$\overline{\tilde{\psi}}^{w} = \arctan \left( \frac{\tilde{\zeta}^{w} \circ \sin \left( \tilde{\psi}^{w} \right) + \tilde{\zeta}^{w} \circ \sin \left( \tilde{\tilde{\psi}}^{w} \right)}{\tilde{\zeta}^{w' \circ} \cos \left( \tilde{\psi}^{w} \right) + \tilde{\zeta}^{w} \circ \cos \left( \tilde{\tilde{\psi}}^{w} \right)} \right) \quad \forall m \in \{1, \ldots, M\} \tag{35}$$

where the symbol represents the *Hadamard* product, and **ζ** *m* and **ζ ^***m* are two *P*-dimensional vectors. In particular, **ζ** *m* = *ζ*<sup>1</sup> *m* , ⋯, *ζ<sup>P</sup> m* T = *ζ*<sup>1</sup> (**ψ˜** *m* ), ⋯, *ζ<sup>P</sup>* (**ψ˜** *m* ) T is expressed in terms of the *spatial coherence* relevant to the *original* multilook interferogram **ψ˜** *m* . Similarly, **ζ ^***m* = *ζ* ^ 1 *m* , ⋯, *ζ* ^ *P m* T = *ζ*<sup>1</sup> (**ψ˜ ^** *m* ), ⋯, *ζ<sup>P</sup>* (**ψ˜ ^** *m* ) T is expressed in terms of the *spatial coherence* relevant to the *reconstructed* multilook interferogram **ψ˜ ^** *m* . The block diagram of the MCh-NF algorithm is depicted in Figure 5. A pertinent pseudo-code for computing the filtered interferometric data stack is also presented (Figure 6).

Note that the exploitation of "conventional" small baseline multilook interferograms is the distinctive characteristic of MCh-NF approach with respect to other previous solutions, such as the *SqueeSAR* [44] and *Phase Linking* [62] methods and other recently proposed multitem‐ poral-filtering techniques [60, 61] based on constraining the analysis to distributed scatter‐ ers [29], which are identified through a pixel-by-pixel selection procedure performed at the full resolution complex SAR image spatial grid. Such a selection permits to rely on the distributed scattering hypothesis, under which the probability density function (pdf) of the complex-valued SAR image may be regarded as being a zero-mean multivariate circular normal distribution, and an appropriate maximum likelihood (ML) estimation step of the filtered phase values associated to each SAR acquisition is implemented. On the contrary, the presented MCh-NF approach focuses on conventional multilook interferograms generat‐ ed without any *a priori* pixel selection step. Accordingly, in this case, it is not possible to rely on the validity of the above-mentioned distributed scattering hypothesis. Therefore, both the phase linking [62] and the phase triangulation of the SqueeSAR [44] algorithms require a preliminary identification of the statistically homogeneous pixels (SHPs) on the full-resolu‐ tion range-azimuth grid. In particular, in ref. [44], the selection strategy of these pixels is based on the application of the *Kolmogorov–Smirnov* test to carefully select a homogeneous statistical population. Clearly, this requires working at the full resolution spatial scale and implies the analysis of the amplitude values of the complex SAR image pixels. A more detailed comparison among the presented MCh-NF method and the ones provided in ref. [44] can be found in ref. [43].

**Figure 5.** MCh-NF block diagram

Figure 5. MCh‐NF block diagram

#### **4. Experimental results 4. Experimental results**

We present in this section some results we obtained by processing a data set consisting of 39 SAR images (Track 308, Frame 2754), collected by the ENVISAT sensor between December 2002 and August 2010 over the Abruzzi region (Italy). The test‐site area includes the city of L'Aquila and its surroundings, which were struck on 6 April 2009, by an Mw 6.3 earthquake that caused more than three hundred fatalities, thousands of evacuees, as well as severe industrial and residential building damages. Starting from the available SAR images, we retrieved a stack of 338 small baseline differential SAR interferograms with a maximum We present in this section some results we obtained by processing a data set consisting of 39 SAR images (Track 308, Frame 2754), collected by the ENVISAT sensor between December 2002 and August 2010 over the Abruzzi region (Italy). The test-site area includes the city of L'Aquila and its surroundings, which were struck on 6 April 2009, by an Mw 6.3 earth‐ quake that caused more than three hundred fatalities, thousands of evacuees, as well as severe industrial and residential building damages. Starting from the available SAR images, we retrieved a stack of 338 small baseline differential SAR interferograms with a maximum perpendicular baseline of 400 m and a maximum time span of 2000 days [43]. The interfero‐ grams have been computed by performing a complex multilook operation with 4 and 20 looks in the range and azimuth directions, respectively. For the interferogram generation, we used precise satellite orbit information and a three-arcsecond *shuttle radar topography mission*

```
 (SRTM) digital elevation model (DEM) of the region to remove the topographic phase
contributions. Finally, the multilook interferograms have been prefiltered by applying the
Goldstein's filter [42].
         Replace with: 
                Topological Characterization and Advanced Noise-Filtering Techniques for Phase Unwrapping of...
                                                                     http://dx.doi.org/10.5772/61847
                                                                                                    363
```


**Figure 6.** Pseudo code of the MCh‐NF algorithm.

To investigate the performance of the presented noise-filtering approach, we applied the nonlinear minimization procedure in Eq. (32) to the stack of the generated (original) multilook small baseline interferograms. As a result, we retrieved a new set of noise-filtered interfero‐ grams that are characterized by generally improved coherence values. This is clearly visible in Figure 7a–f, where we compare three unfiltered (left side) interferograms with the corre‐ sponding (right side) noise-filtered interferograms. It is evident how the fringes due to the earthquake are well recovered. Such interferometric data stacks can then be used for the generation of surface deformation time series using the small baseline subset (SBAS) [6]

**Figure 7.** Comparison between the original (left column) and noise-filtered (right column) multilook interferograms relevant to the area of the Abruzzi region (Italy). (a–c) October 30, 2005, to November 8, 2009, August 21, 2005, to June 6, 2010, and August 1, 2004, to April 12, 2009, interferograms, characterized by perpendicular baseline values of 192, 145, and 395 m, respectively. (d–f) Noise-filtered multilook interferograms corresponding to the ones in panels a–c, re‐ spectively.

processing chain, whose parallel version (P-SBAS) has been proposed in refs. [13, 14, 35, 36]. This is matter for the analysis presented in the next subsection.

Topological Characterization and Advanced Noise-Filtering Techniques for Phase Unwrapping of... http://dx.doi.org/10.5772/61847 365

**4.1. The use of MCh-NF algorithm within MCh-PhU framework**

We present in this subsection how the MCh-NF algorithm can be efficiently used within the SBAS processing chain, where phase unwrapping procedures are implemented through the MCh-PhU technique known as *extended minimum cost flow* (EMCF), also discussed in the first part of the chapter. Figure 8 shows the diagram block of the advanced EMCF-based SBAS processing chain [26, 35, 43], which integrates the conventional SBAS codes, exploiting the EMCF MCh-PhU procedure to perform phase unwrapping operations, with the presented MCh-NF noise-filtering technique. In addition, an effective procedure for the selection of timeredundant interferograms is also included; interested readers can find additional details in ref. [43]. To provide an example of the potential of the advanced processing chain incorporating both the discussed MCh-PhU and MC-NF techniques, we here focus on the area of *Yellow‐ stone* caldera, representing one of the largest and most active volcanic systems in the world. We analyze the temporal evolution of the surface deformation occurring in this area by applying the implemented EMCF-SBAS processing chain to a set of 22 ENVISAT images (Track 41, Frame 2709), acquired from May 2005 to September 2010, from which we have retrieved a corresponding set of 122 small baseline interferograms [43]. As in the previous case study (relevant to Abruzzi area), the prescribed maximum values of 400 m and 2000 days for perpendicular and temporal baseline, respectively, have also been considered. The retrieved mean deformation velocity map depicted in Figure 9, which has been obtained by applying the processing chain (MCh-NF + EMCF-SBAS) of Figure 8, allows us to recognize the complex deformation scenario affecting the Yellowstone Caldera region and its surroundings, where uplift effects and broad subsidence patterns characterize the detected deformation field. In addition, the deformation time series relevant to four pixels, whose locations are identified by the black stars labeled as P1, P2, P3, and P4, are also illustrated in Figure 9.

## **5. Conclusion**

Within the context of SAR interferometry, two main issues related to the multichannel phase unwrapping and noise filtering for interferometric data stacks processing have been ad‐ dressed. First, a rigorous gradient-based formulation for the multichannel phase unwrapping (MCh-PhU) problem has been systematically established, thus providing a topological

**Figure 9.** (a) Mean deformation velocity map (in color) of Yellowstone Caldera, computed in coherent pixels only and superimposed on the SAR amplitude image (gray-scale representation) of the zone, retrieved by applying the ad‐ vanced EMCF-SBAS processing chain. The black square marks the location of the reference SAR pixel. (b–e) Deforma‐ tion time series relevant to the four pixels identified via black stars in Fig. 8(a).

characterization of the problem within the purview of the theoretical foundation of the *discrete calculus*. Then the innovative MCh-NF procedure for the noise filtering of time-redundant multichannel multilook interferograms has been properly presented within the considered topological framework, by adopting a consistent formalism. Finally, some experimental results obtained with real data have been shown, thus demonstrating the effectiveness of our approaches and their relevance for geospatial phenomena analysis and understanding.

# **Acknowledgements**

This work was supported by the Italian Ministry of University and Research (MIUR) under the project "Progetto Bandiera RITMARE." We would like to thank the European Space Agency for providing the ENVISAT ASAR data and the University of Delft, Delft, The Netherlands, for the related precise orbits. We would also like to thank Italian Space Agency (ASI), which has provided us the Cosmo-SkyMed SAR images under the framework of the European Union's Seventh Program for research, technological development, and demon‐ stration MED-SUV project (grant no. 308665). Finally, the authors thank the Japanese Space Agency (JAXA), which has provided the used ALOS-1 data through the project entitled "Advanced Interferometric SAR Techniques for Earth Observation at L-band" (ID project 1149) in the framework of "The 4-th ALOS Research Announcement for ALOS-2" call.

## **Author details**

Pasquale Imperatore\* and Antonio Pepe

\*Address all correspondence to: imperatore.p@irea.cnr.it

Istituto per il Rilevamento Elettromagnetico dell'Ambiente (IREA), National Research Council (CNR) of Italy, Napoli, Italy

# **References**


metric phase images," IEEE Trans. Geosci. Remote Sens., vol. 36, no. 5, pp. 1456– 1465, Sep. 1998.

